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Abstract. Many physical systems can be adequately modelled using a second order approximation. 
The problem of plant identification reduces to the problem of estimating the position of a single pair 
of complex conjugate poles. 

One approach to the problem is to apply the method of least squares to the time domain data. This 
type of computation is best carried out in "batch" mode and applies to an entire data set. Another 
approach would be to design an adaptive filter and to use autoregressive, AR, techniques. This would 
be well suited to continuous real-time data and could track slow changes on the underlying plant. 

I this paper we present a very fast but approximate technique for the estimation of the position 
of a single pair of complex conjugate poles, using the Hilbert transform to reconstruct the analytic 
signal. 



INTRODUCTION 

In the theory of control, it is most common for physical systems to be mathematically 
modelled using coupled linear systems of ordinary differential equations [JIp. When these 
equations are transformed using integral transforms, such as those of Laplace or Fourier, 
then the physical systems are modelled using finite rational polynomials in an auxiliary 
variable, s = jco: 
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The zeros of the polynomial, Q(s), are called "poles" and correspond to responses that 
have finite output for zero input. These are called "modes." It is very common for one 
mode to dominate the response of the whole system. It is also common for this mode 
to be of a damped oscillatory type, corresponding to a single pair of complex conjugate 
poles. This can occur whenever the potential energy function of the system has a local 
minimum [0 In this case, we can approximate a large complicated system, with many 
poles and zeros, by a simple second-order system with a single pair of complex conjugate 
poles. This is called the second-order approximation. Many mechanical or electrical 
systems can be realistically modelled using the second-order approximation. We can 
write: 
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If we wish to model the behaviour of a real physical system, using an approximate 
second-order model, then it is necessary for us to estimate the position of the pole pair. 



This could be done in the frequency domain, by exciting the system with a sinusoidal 
source and then measuring the magnitude and phase of the response at different fre- 
quencies, but this is often not practical. There are situations when the only practical 
sources are step functions, l(t) or impulses 8(t). We can excite the system with steps or 
impulses and then sample the response in the time domain. The impulse response of a 
second-order system will generally be of the form: 

y(t) =Ae~ at cos (co d t) +Be~ at sin (co d t) (3) 

where a> d z = co 2 - a 2 . 

The problem of plant identification then becomes equivalent to asking: How do we 
estimate the position of the pair of complex-conjugate poles if the only data at our 
disposal is a set of time-domain samples of the response of the system to steps or 
impulses. 

As a possible illustration, we could imagine that we strike a bell with a hammer and 
then record the sound as it gradually decays. We want to estimate the damped frequency 
of oscillation, 00 d , and the damping coefficient, a, using only the data from our sound 
recording. 

If we knew something about the distribution of the errors of measurement then we 
could apply the method of maximum likelihood to estimate the parameters , co d and a. 
If the errors were known to be the result of a very large number of uncorrelated random 
effects then we could apply the Central Limit Theorem and we could assume that the 
errors had a Gaussian distribution. The problem of plant identification would reduce to 
a non-linear least-squares estimation problem 0]. The difficulty with this approach is 
that the resulting equations would be non-linear and would have to be solved iteratively, 
using a numerical method such as gradient descent. A further weakness of this approach 
is that it would be an exact solution to an approximation of the real problem. It would 
be far more reasonable to have a quick but approximate solution to the approximate, 
second-order, problem. This would tell us most of what we need to know without having 
to waste a lot of effort. 

In this paper, we present a fast, but approximate, algorithm for the estimation of 
the position of a complex conjugate pair of poles on the s plane. We use a discrete 
approximation to the Hilbert Transform , to reconstruct the complex analytic signal 
form the sampled real time signal. The analytic signal has a complex-exponential form. 
We apply very simple statistical techniques to the analytic signal in order to obtain the 
required parameters. This approach is an alternative to the more conventional, non-linear 
least squares or Autoregressive, AR, approaches to the problem. 

The problem of parameter identification, for a freely vibrating system has been studied 
by Feldman [Q] who used the Hilbert Transform to provide information about instanta- 
neous amplitude and phase of a signal. The method that we present here is more simple, 
and limited, than the approach used by Feldman. 



THE ANALYTIC SIGNAL OF THE SECOND ORDER RESPONSE 



The response described in Equation [3| is equivalent to 

y(t) = Ce~ at cos ((D d t - (j)) (4) 

where C = VA 2 +B 2 , cos(» = A/VA 2 + B 2 and sin(» = B/VA 2 + B 2 . This type of 
function will apply whenever the input to the system is zero. If the input is a finite sum 
of step and impulse function then the input will be zero for most of the time. There will 
be abrupt changes in C and but the parameters, a and co d will be constant as long as 
the structure of the plant is maintained. 

The immediate aim is to reconstruct the analytic signal. The Hilbert transform ^ [7|] 
is a standard technique for achieving this. The Hilbert transform of u(t) is defined as: 

v{s)= ir^i dt . (5) 
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It has the important property that it defines the relationship between the real and imagi- 
nary parts of a complex analytic function. If we have an analytic function: 

<t>(t) = u{t)+jv(t) (6) 

with real and imaginary parts u(t) and v(t) then the relation for v(s) is given by Equation 
|5[ Bedrosian's theorem tells us that the Hilbert transform of a(t) cos(fttf) is a(t) sin(ft>f ). 
If we apply Bedrosian's [jfjtheorem and the shifting property to Equation |]then we find 
that the analytic signal is: 

y 2 (t) = Ce~ m cos (o d t - 0) +jCe~ at sin (ft)/- 0) (7) 
= Ce- j ^e^ a+ja d)' . (8) 

This analytic signal is a pure exponential function and can essentially be "unwrapped" 
using the log() function. We can write: 

log (y 2 (t)) = log(C) - ;0 + (-a + j(O d ) t . (9) 

If we sample this analytic signal at intervals of T s then we can numerically calculate the 
slope of Equation |9| to get: 

_ a+]aj= ^M±J^iMA\ (10) 

This allows us to directly estimate the parameters, a and C0 d . We note that the use of 
the difference operation has removed all reference to C and to 0. This means that the 
method is not sensitive to the initial conditions that apply immediately after the shocks 
that occur then the impulses and steps are fed into the system. Our only requirement is 
that the input to the system is zero for most of the time. 



The Matlab code required to implement this algorithm is very short and simple: 

% reconstruct the analytic signal, y2 
y2 = hilbert ( y ) ; 

% take the natural logarithm, log(y2) 
L = log (y2 ) ; 

% unwrap the phase of the log(y2) 

L = real (L) + j *unwrap ( imag (L) ) ; 

% estimate the differences of the log of y2 

D = cliff (L) ; 



A SIMPLE STATISTICAL TECHNIQUE 

Equation WU suggests that we should be able to precisely estimate the required parame- 



ters. There are a few practical problems with the direct application of Equation |T0|: 

• The numerical calculation of the Hilbert transform relies on the Fast Fourier 
Transform and there are limitations imposed by the finite number of samples. These 
includes the "Gibbs effect," due to the finite length of the data set. 

• Real samples from a physical process will be subject to noise and errors of 
measurement. The differencing operation tends to magnify the effect of noise. 

• There would be a number of outliers caused by the "shocks" of the steps or 
impulses. Some measurements will not be reliable. 

The authors have found that the median is a very robust measure of location and is 
less subject to the influence of the outlying values than the arithmetic mean. The Matlab 
code for this is very simple: 

% calculate the real and imaginary parts of the differences 
d_mag = real (D) ; 
d_phase = imag(D) ; 

% calculate the median rates of change 
mid_re_slope = median (d_mag) ; 
mid_im_slope = median (d_phase) ; 



SOME RESULTS 

A second order system was simulated, using known parameters, and the parameters were 
then estimated using the new algorithm. The reconstructed analytic signal is shown in 
Figure [I]. 

This same data can be represented in three dimensions. This is shown in Figure [Z[ 
The logarithmic slopes were estimated and the parameters were calculated. 
The quality factor of the simulated system was Q — 10 and so the theoretical value 
of the damping constant was a = 0.05. The estimated value was a £St = 0.0472. The 
characteristic frequency of the simulated plant was co d = 1.0 and the estimate was 
m , = 0.9948. 
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FIGURE 1. The original sampled signal is shown on the left. The signal on the right is an estimate 
of the complex part of the analytic signal which was reconstructed using the Hilbert Transform. There is 
significant error in the reconstructed signal near the sample boundaries, due to the Gibbs effect. This is 
the result of fact that Matlab uses the FFT to calculate the Hilbert transform. 



There was a —5.5% error in the estimate of the damping constant, a and there was a 
—0.4% error in the estimate of the damped natural angular frequency of oscillation, co d . 



SUMMARY AND LIMITATIONS 

The method does work as an approximation and could estimate a to within —5.5% and 
(O d to within —0.4%. The error in the estimate of the damping constant, a and there was 
a —0.4% error in the estimate of the damped natural angular frequency of oscillation, 
(O d . There are some notable problems with the method: 

• It is sensitive to relative time scale of the constants in question and the width of the 
sampling window. 

• It suffers from boundary effects due to the Gibbs' phenomenon. 

• The method is also numerically unstable for large data sets with large values of 
time, t, since exp(-at) can cause Matlab to underflow. 
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FIGURE 2. This is a 3D plot of the reconstructed analytic signal. The XY plane and all planes parallel 
to it represent the complex field that contains the analytic signal. The vertical, or Z, axis represents time. 
The general appearance of damped oscillation is unmistakable. 



Our simulations suggest that this is a fast but approximate technique for the estimation 
of the position of a single pair of complex conjugate poles but it is sensitive to a number 
of factors which may limit its practical application. 
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